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Tice ABSTRACT 


In this study, the flow characteristics of air arcund a 
Gateular cylinder were determined by two means: by 
analytical solution of the Navier-Stokes eguations, and 
actual measurement of the flow itself. A DISA LD 5586 laser 
doppler anemometer system was employed to measure normal (v) 
and parallel (u) velocity components of a laminar flow 
regime about a three-sixteenths inch diameter cylinder 
mounted in a 32 x 45 inch wind tunnel (Reynolds number 


approximately 40). 


A finite element grid was constructed using triangular 
elements which encompassed the upper half-cylinder. Flow 
measurements using the laser velocimeter were confined to 


this region cf interest. 


Theory Was compared against experimental results, and 
the feasibility of the utilization of this particular type 


of laser anemcmeter was evaluated. 
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ens DESCRIPTION OF LASER ANEMOMETER 


ae SS = we sli eee 








A basic laser anemometer is comprised of four primary 


elements: 


(1) a coherent light source 


(2) optics 


(3) Signal converter (photomultiplier tube) 


(4) signal processor/display 


The principle of the laser anemometer is to create a 
grid pattern and measure the speed with which a particle in 
a fluid free stream passes through the grid. The following 


assumptions were taken into consideration: 


(1) particles injected into the free stream are 
small enough to faithfully follow strean 
lines of the flow, thus by measuring particle 
characteristics, we are actually measuring 
the medium's characteristics. This is easily 
achieved by use of particles which are on the 


order of a few microns in diameter. 


(2) particle generation results in a consistent 
particle size, or at least a majority of the 
particles are within a specified range of 
sizes. For example, DOP (di-2-ethylhexyl 
phthalate) particle generators produce 
particles with a mean diameter of Os 5 
microns and burning turbine oil precduces 
particles on the order of 4-6 microns (from 


electron microscope). 





A laser is incorporated because of its coherency and 
power. The aforementioned grid is in reality the result of 
two coherent beams interfering with themselves (Fig. 1). 
The optics produce the beam pair by conventional splitting 
and focusing. Because the laser beam intensity is of a 
Gaussian distribution the grid created by constructive and 
destructive interference is not uniform over the entire 
focal volume. The fringes created on the outer porticns of 
the horizontal axis are not perfectly formed because of 
Slightly different intensities in that area. The main 
portion of the interference grid is uniform and, most 


importantly, the fringe spacing Oe 1s determined only by 
the laser wavelength (lambda) and the half angle of the bean 


pair convergence (theta/2). 
e = lambda/ (2*sin (Theta/2) ) 


The number of fringes eee may also be comruted, 
provided the beam focusing convergence angle (delta theta) 


is known ,by the relation: 


te cone 474). 7 (PTsDelta Theta) 


The number cf fringes is important only to ensure that the 
signal processing apparatus will receive enough information 
to . make possible further data reduction. If an insufficient 
number of fringes exists, not enough signal will be 


available for velocity determination. 


Solid particles with a mean diameter which is less than 


half the fringe spacing are introduced into the free stream 
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of the fluid to he measured. When the particle is small 
ccmpared to the fringe spacing the scattered light exhibits 
essentially 100% modulation, while particles equal to the 
fringe spacing produce negligible modulation{Ref. 2]. The 
particles may be produced by atomization of water, dop-air 
reaction, burning oil, smoke bombs, or may be kought 
commercially, as e.g. polystyrene spheres and some _ solid 
chemicals such as zinc oxide. Choice of the particle to be 
used is normally dictated by the fringe spacing in the focal 
volume, and the nature of the medium to be measured. 
Regardless of choice, the particle must be small enough to 
follow the streamlines as closely as possible, because the 
laser anemometer actually measures the particle velocity and 
unless the particle's motion is a true representation of the 


fluid velccity, pertinent data may not be obtained. 


Inside the measuring (or focal) volume, bands of 
darkness and brightness exist due to the interference cf the 
two beams. Particles entering the focal volume will 
alternately scatter light as they pass through regions of 
constructive interference ( typical focal volume dimensions 
are 300 microns wide by 4 millemeters long and contain 
approximately 64 fringes). A photomultiplier tube located 
on the optical axis and focused into the measuring volume 
receives this scattered light and converts what appears to 
be a modulated light signal into a series of voltage pulses. 
Each voltage pulse correlates to a particle crossing one 
fringe. Fecause the spacing of the fringes is known, the 
frequency of the voltage pulse train may be equated to the 
velocity of the particle by the following relationship: 


£=V ae = 2*V *sin(Theta/2) /lambda 
oO oO 


where oy is the fringe spacing and uP is the velocity of the 








particle. Now the processor is left with a fairly simple 
task of ccunting the frequency and multiplying it by a 
calibration factor, which is, in essence, a fringe spacing 


parameter to arrive at a velocity. 


As long as we are assured that only one particle is in 
the measuring volume at one time, the apparatus will operate 
correctly, but this is an unrealistic assumption. [In order 
to achieve such a situation, the seeding density of the 
particles would be exceptionally sparse and the resulting 
data rate would be- so slow as to preclude timely data 
reduction. To overcome this obstacle, a validation circuit 
is incorporated into the counter/processor. When a particle 
crosses a fringe spacing, two counters are triggered. One 
measures the time duration for five fringe crossings and 
holds this information. The second counter continues to 
measure the time for a total of eight fringe crossings. The 
time of the first counter is then compared with that of the 
second counter, and if the first time is indeed five-cights 
of the second time (within specified limits, i.e., comp. 
accuracy) a validation gate is opened and the resultant 
frequency is then stored into an ensemble cell. The ensemble 
cell is adjustable to hoid various numbers of counts and 
average them for the final velocity display. The final 
displayed velocity then is that component of fluid velocity 
which is perpendicular to the optical axis and contained in 


the same plane as the beam pair (Fig. 2). 


The anemometer May be operated in three _ modes: 
reference beam, forward scatter, and back scatter. The last 
two modes are most frequently employed and are described 
here: 

Forward scatter mode: in this mode, the photomultiplier 
tube is located on the optical axis and focused in a 
direction directly opposed to the direction of the emerging 


beam pairs. This mode offers the strongest attainable 
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Signal. 

Back scatter mode: again the photomultiplier tube is 
mounted on the optical axis, but now it is aligned in the 
direction of the emerging beam pairs. This mode affords the 
ability tc measure fluid velocities in confined spaces or in 
Situations where the dectector cannot be mounted in fcrward 
scatter, such as in measuring velocities in a compressor. 
Additionally, operation in this mode requires a higher 
powered laser due to the reduced intensity of back 
scattered light (Fig. 3). 
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III ANALYTICAL DERIVATION 

The finite element method is a numerical analysis 
technique fcr obtaining approximate solutions to a wide 
Variety of enginering problems. [It offers a means to solve 
complex continuum problems by sub-dividing them into simpler 
interrelated problems. The nature of this engineering 
problem does not readily lend itself to a closed fora 
solution. Therefore a numerical solution is dictated. 


The method of finite elements is well employed where 
there are irregular geometries or unuSual specifications of 
boundary conditions. The method gives a piecewice 
approximation to the governing equations. Figure 4 shows 
the finite element grid employed in this flow problem. 
Approximating functions (interpolation functions) are 
defined in terms of the values of the field variables at the 


nodes. = 


For the finite element representation, the nodal 
values of the field variables hecome the new unknowns and, 
once found, the interpolation functions define the field 
variable throughout the field of elements. The degree of the 
interpolaticn polynomials is governed by the number of nodes 
per element and dictated by the order of the governing 


equations. 


The basic principle of the finite element method is to 
divide the solution domain into a finite number of 
subdomains (elements) which are connected only at node 
points. Interpolation polynomials are created for each 
element anda collecton of interpolation functions over the 
whole soluticn domain provide a piecewise approximation to 
the field variable. 


Introduction of generalizations facilitate direct 


derivation of finite element equations from governing 


differentiai equations. The method of weighted residuals 
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(Gaierkin's Method) is a technique for obtaining approximate 
solutions to non-linear partial differential equations. We 
first assume the general function behavior of the dependent 
field variable in some way so aS tO approximately satisfy 
the given differential equations and boundary conditions. 
After applying this approximation to the original equations, 
scme error will result, which we desire to minimize or drive 


to zero over the entire solution domain. 


The basic set of differential equations governing the 
“problem from a velocity and pressure formulation approach 
are (after linearization by approximating non-linear viscous 


terms u,v): 
n n 


ae Ou 
eonttnuity) Bx ~ 5, = O 


Lin ~*~ — -+ gP e 
ncmentum | Ox Va ay e vie Vy Ce 


u, v, and EF are approximate solutions to the flow problen 
n n n 


velol 
and may be visualized as the hn iteration of the procedure 


OL SOlUtLOn . These equations used in formulation of the 
problem are by no means unique. Stream functions and stream 
function with vorticity formulations have also been employed 


to establish finite element solution procedures. 


Inasmuch as the velocity terms are of second order, a 


second order interpolation polynomial must be used to 
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Calculate these unknowns. The pressure terms are of first 


order and hence reguire only first order interpolation. 


The basic element chosen for this solution was the 
triangle and, in a localized scheme, is visualized as 


follows: 


Nodes 2, 4 and 6 lie equidistant between their respective 
corner nodes. 


Nodes 1, 3, and 5 are referred to as corner nodes and 
are the basis for the linear interpolation polynomials 
derived fcr the pressure terms. For any triangular element, 
the corner nodes are described by their Cartesian 
coordinates Xx , y..- 

mt 1 
Consider the following local element: 


rY_) 





The coordinates of node one are xs geo node two xX, y $3 
etc. and the interior point (x,y) may be described as a 


linear combination of the three corner nodes: 
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Zee eri xs tL x 
1° Si eva 3 3 

=L y +L y +L 
Y ey ae ae 


where oF is the natural coordinate (interpolation 
polynomial, here linear) or weighting tunction which relates 
the coordinates of the corner nodes to any interior node. 
Le has a value of one at the point x OY, and is zero at the 
other two nodes of the element. 


The L 's expressed as linear functions of the Cartesian 
a 


coordinates then are 


\ 
ee) - ez a) tc Y) 
i eae % 
ae ty) ZA (ape Xtc) 


pee 
Lj (zY)= ZR (a+b xtc y) 


where 2A = PY 


which is equal to the area of the triangular element and 
a =x =X b =y - and c =y - the determinants of 
1 *2t3 *3%2! 147273, *% 1%, 72° 
the minors corresponding coefficients aus Die orc and the 
other coefficients are similarly computed. 
For the velocity terms, second order interpolation 
polynomials are indicated, while the pressure terms remain 


first order. Hence the introduction of more nodes (6) will 


result ina new set of weighting functions which may also be 
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expressed in terms of the natural coordinate interpolation 


pclynomials: 

N (=2L"-L N =4L L 
1 2 1. ce 
a 

N 37 <r N =4L L 
2 4 3 
Ze 

N 57h N =4L L 
3 6 i 3 

P 

N =L 

me 4, 
6 e 6 e 3 ?P 

u= > Nu. v= DNV. P= DN P. (xxxx) 
ea 1 ts ae ’ A 


Application of Galerkin's Method to the governing 
equations results in the following intergral expressions: 
Su Sue eee IGS CO 
SNe ay * V, vu) 
( if ox Bey ” P ox 


3 HO ee lene bs 
G Ni (un 5S ew al Vv?v) do. O 


(P) 
OV , 
here 2G\e implies "over the area of the triangle". 


Then we may make the substitution of the field variable 


in terms of x and y into terms of L.. 


ai 
ergy: given u (X,Y) => u(L) 
i 

then ou. eb IN, JL 

Ox OL; ax 
where Li = bi 

ax = 2A ay. w aNi abi 
and a similar procedure for ay > Uo ali dy 
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where 


application of the above chain rule allows integration by 


Parts: 
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The evaluated integral expressions tnen form the influence 
Matrix as fcllows: 
i l 
K :0;,K 
ae ee 


O,K'K 
Ly; 3 


~-p 4-3 


KS 1K 10 
L =) 3X LJ 


Specification of the field variables or right hand sides of 
the matrix eguation then provides a defined problem: 
(r4- [x] = fo] 
Additionally, evaluation of the integral expressions yields 
boundary terms which must be computed for nodes which lie on 


the boundary of the control surface: 


* 46 
R = ( N X ds R = ( N Y ds 
1 i 2 i 
* A A 4 @ain pe 
where X = uyuen-h, Pe and Y = 4TV T= ny P 


It should be noted that along vertical boundaries there is 
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no outward pointing normal n and Slmilarly along horizontal 
ag 


e A e 
boundaries no n exists, and furthermore these terms are not 
x 


computed where u and v are specified. Additionally, it was 
assumed that along the left-and right-hand side and upper 
boundaries ®%/ax and “fay were identically equal to zero 
because these boundaries were placed at five’ or nore 
characteristic lengths (diameters) from the cylinder and 
considered to be outside the realm of influence of the flow 


about the cylinder. 
The terms comprising Ke K and KS are entered element 


by element into the main program in steps FLU02080 through 


QZ00170. Then each element's matrix is compiled into the 
field matrix (comprised of all elemental matricies) and is 
solved iteratively to output the horizontal velocity 
component u, the vertical component v, and the pressure P at 


every node in the field. 
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Iv. | DESCRIPTION OF EXPERI 
A right-angled, three-sixteenths inch cylinder made of 
stainless steel was mounted in an AEROLAB series 90 32x45 
inch low speed wind tunnel through the floor of the test 
SCGetemeitig. 59). The portion of the rod protruding down 
through the test section was mounted on a 
two-degree-of-freedom, micrometer controlled bed which 
allowed the rod to be maneuvered in the vertical plane 
parallel to the flow direction. This scheme enabled the 
laser velocimeter and photomultiplier tube to be nounted, 
aligned and left undisturbed for the entire period during 
which data were taken. The bed waS positioned to read 
Locations in the control volume which directly coincided 
With the location of the nodes in the finite element mesh 
(i age ©} ee 


In the interest of precluding as many opportunities for 
error aS possible, all u velocity components were first read 
and then ccmponents in a forty-five degree position were 
taken thus necessitating only one re-orientation of the bean 


pairs. 


The plexiglass window prefacing the optics transducer 
of the test section had one-quarter inch holes drilled into 
1t to allow the beam pairs to pass unimpeded, thereby 
Minimizing distortion of the focal volume fringe pattern. 
The photomultiplier tube was likewise accommodated through 


the far window. 


The seeding mechanism employed Zinc-hexaclorcthane 
based (non-toxic) white smoke bombs. The nature or the 
closed return wind tunnel permitted continuous measurement 
Without interruption for re-seeding of the flow. By means of 
electron microscopy, the average size of the smoke particles 


was determined to be on the order of two microns (Fig. 7). 








Cylinder mounted in wind tunnel and LDV setup. 


Figure 5 
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Zinc-hexachlorathane particles 
(X 5000 lor= 2.0 um) 


Figure 7 
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The final velocity component measured at each node was 
the result of averaging three sucessive readings from the 
velocimeter processor. Each of these readings represented 
the average of 4096 readings which had been validated to 
within 1.5% of the exact five-eighths validation count 


criteria. 


The flow in the wind tunnel was stabilized at five 
inches per second to maintain a Reynolds number of about 40, 
which is in the realm of creeping flow. The experimental 
setup was designed to measure the steady rlow about the 
cylinder, hence the relatively low Reynolds number, and no 
effort was made to obtain unsteady flow measurements. Three 
separate runs were conducted, two over the forward half of 
the contrel volume and one over the rear half. Because the 
beams of the laser anemometer were approximately 
one-sixteenth inch in diameter, several node position 
measurements were excluded due to inaccessibility. The 
majority of these were in close proximity to the cylinder. 
Mm order to obtain both an x- and a y-velocity component of 
the flow at every node which was accessible, two 
measurements were made, but because of the very small 
vertical ccrronent present in most of the flow, direct 
vertical beam pair measurement was unsuitable. Instead, 
after obtaining the horizontal component for each node in a 
run, a second series of measurements was made with the bean 
pairs inclined forty-five degrees from the vertical. In this 
fashion the y-component was determined by the following 


relationship: 








It can ke seen here that the unknown resultant (R) aas 
a projection onto each axis representing the orientation of 
the beam pairs, thus the resultant is the vector which 
emanates from the origin and terminates at the intersection 
of the perpendiculars constructed from the uU and v(at 45 


degrees) vectcrs. 
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OMPUTER PROGRAM 





The program found in the appendix has several ccmment 
cards inserted for clarification. However, a broader 
presentation is offered here. The program may be subdivided 
into four main sections: 

Section 1: dimensioning and zeroing of 
parameters and arrays. (FLUO0020 through 
FLU00690) 

Section 2: input of node coordinates, system 
topology and known values of uu, v, P, and 
right hand sides, and output of the same. 
(FLU00700 through FLUO2090) 

Section oe computation of interpolation 
polynomials on an elemental level and 
arrangement into Local element arrays. 
(FLU02100 through GO00320) 

section 4; loading of the global influence 
matrix (TM) and solution of (TM) (T)=(Q) by 
Gauss-Jordan partial pivoting and back 
Substitution, with an iterative routine 
incorporated until convergence is attained, 
and associated print out. (FLU0O4470 through 
FLU04880) 


The computer routine systematically tested the solution 
of field variables with the previous solution until all 
terms in the T and T1 arrays agreed within 0.001 (ESPILA). 
If the convergence test failed, the process was repeated 


from section 3. 


Initially, the program was run on an IBM 360 in FORTRAN 
CLG mode, then an ensuing run was completed on the FORTRAN 
HCLG compiler and a time reduction of two-thirds was 
realized. The program, as given in the appendix, required 


approximately 375K bytes of core to compile and 675K bytes 
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to execute and output. Substantial space savings may be 
realized by utilizing banded matrix storage methods vice the 
method employed here, but due to time limitations this 


avenue was not fully explored. 


The program attempts to solve simultaneously for all 
the field variables. However, introduction of zeros on the 


main diagonal of the TM influence matrix resulted in a 


Singular matrix. Another approach was explored where the 
equation of continuity 

Ol 4. OV = 

ax * ay 0 


was replaced fry 


-ty'P - ie - uy | 


and the resulting terms were added to the TM matrix and 
their respective right hand sides. The only field variables 
affected by this form of the continuity equation were the 
pressure variables, and depending on whether the lacal 
corner node involved in the local pressure equation was on 
the boundary or interior, different right hand sides had to 
be evaluated. This is the purpose of steps FLUO4460 through 
GO00320. The Simultaneous solution technique for all field 
variables applied to this new set of equaticns did not 


converge. 


A third, unimplemented technique may be considered. 


Given the original matrix equations: 


(1) KU + 0 +K P =R 

a 2 1 

(2) 0 +K V +K P =R 

1 3 2 

(3) K U +K V + O =0 
2 3 


Solution for u and v will result in 


ane =-1 
(1) U=-K KP +K_ R 
1 2 1 1 


Zo 





=-1 -—1 
(2) V=-K KP +#K- R 
1 s 1 2 


and multiplication by K and K yields 
ee —1 
K O=-K K K P+ K K R 
2 aael 2 2 
-1 -1 
K V=-K K Keep ft OK K R 
3 3 3 
Summation of these equations allows substitution for the 


third matrix equation: 


—1 -1 -1 -1 
[k K Kh +k K K ] P= KK Fee ater ae R 
2a 2 3 1 3 


The matrix operating on P may be calculated cn an 


algebraic level and entered, in terms of a, b. and c., into 
1 et 1 


the program in order to solve for the pressures at every 
corner node. The calculated P's then would be entered into 
the field variable array and the partial pivoting routine, 
which involves only those terms where the field variable is 
not specified, would operate primarily -on the first two 
matrix equations. In previous computer trials where all 
cecrner nede fressures were specified, although not known to 
be correct, calculations were properly performed. This 
method would allow the actual pressures present in the 
control surface to be entered into the computer routine for 


further calculation of the velocity components. 
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No computer routine output was obtained for comparison 
of results, due to the inability of the program to converge. 
As such, ccmparison of the two methods cannot be offered 
here. The non-convergence of the program was investigated 
from several aspects: validity of the non-linear terms, 
correctness cf the loading of the TM matrix on both the 
elemental and global levels, correctness of input, 
performance cf the partial pivoting subroutines and the 
iteration scheme; all without success. In test cases the 
program initially solved for the unknowns within magnitude 
accuracy for a majority of terms. However, sucessive 
iterations resulted in oscillating solutions, even after 
twenty iterations. Convergence was expected after ten to 
fifteen iterations without under- or overrelaxation factors 
incorporated into the iteration scheme. After the new 
formulaticn of the continuity eguation was incorporated into 
the computer routine, the program execution again resulted 
in a singular matrix when only the left and right hand and 
upper control surface boundaries had their respective 
pressures specified. Underflow occurred in the execution of 
the subroutine SOLVE at step S0OL00260. Subsequent 


th th 
investigation revealed that the 178 and 179 crows of the 


reduced matrix equations were the cause of the underflow. 


The data which were taken for runs number two and three 
are displayed in the appendix and are plotted on the finite 
element grid in Figure 8. Without the analytical results to 
compare the measured data only a qualitative statement nay 
be made pertaining to the LDV output. It is readily seen 
that the trend of flow directions around the cylinder 
appears to be correct. By and large, the flow in close 
proximity te the cylinder was accelerated over the tcp and 


decelerated pear the stagnation points. It is pointed out 
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here that the plot 1s a composite of two different trials, 
and as such, different free stream velocities were present 
for each run. In Figure 8 the vector representing the flow 
at a node was plotted and the magnitude of the vectcr in 
inches per second was displayed above the vector or by the 
node from which the velocity vector emanated. 
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Velocity data plot 


Figure 8 





A. ERROR ANALYSIS 


Re-positioning of the cylinder to its base location was 
accomplished with relatively low position error. Prior to 
each run the cylinder was aligned uSing the beam pairs as 
reference. In the horizontal mode, the beam pairs would 
strike the leading and trailing portions of the cylinder and 
reflect to positions horizontally opposed on the far side of 
the wind tunnel .wwhen correctly zeroed vertically. The same 
procedure was repeated for the vertical beam pair except 
position was checked by the upper beam reflecting to a 
vertical position on the far wall and the lower beam {at 
reduced intensity) striking the cylinder Support on 
centerline. This procedure was repeated several times and 
the maximum deviation between any aligned positions was 0.2 
millimeters horizontally or vertically. Thus, it was 
determined that the maximum positioning error inherent in 
the system would be the ratio of the smallest coordinate of 
any node measured and the aforementioned 0.2 millimeters or 

022/73. 75=0.5% 


Freestream velocity measurements were likewise tested 
for fluctuations. Typical measurement variations ranged 
Poem O.79129 to 0.142 (run humber three; freestream 
measurements) in the horizontal mode and 0.0947 to 0.100 in 
the forty- five degree mode (velocities in meters/seconds). 
These figures represent a total of 10% variation which may 


be considered as + 5.% from the mean. 
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The experiment, as described, required approximately 
two to four hours for each fun. aside from the initial 
alignment of the components and the cylinder, this was 
perhaps the greatest drawback of the systen. Next, the 
exclusion of regions in close proximity to the cylinder from 
measurement was deemed to be a major factor in the small 
amount of data actually taken. The latter may be improved 
with a different experimental arrangement; particularly, the 
plexiglass section directly prefacing the optics transducer 
Should be completely removed, thus allowing all angles of 
beam pair inclination, and in turn availing amore locations 
to measurement. Another improvement which could oe 
implemented would be unitization of the laser, transducer 
and photomultiplier tube arrangement. This would decrease 
the LDV's sensitivity to vibration and provide a more 
accurate optical axis alignment. The aforementioned 
improvements would still not decrease the time required for 
data collection, but the introduction of a two component 
System [Ref.4] would allow simultaneous measurement using 
orthogonal beam pairs. Furthermore, this type of system 
would not cencern itself with the errors introduced in bean 


pair re-alignments performed with the single beam pair LDV. 


The reproducibility of the experiment was judged to be 
good with the exception of exact duplication of the wind 
tunnel velocity between successive trials. This, of course, 
is primarily a function of the controllability of the tunnel 
motor and fropellor assembly. It is estimated that with 
adequate facilities for incorporation of some of the 
aforementioned improvements, more accurate data may be 
obtained. The LDV displays great potential as a means of 
measuring a wide variety of dynamic fluid regimes accurately 


and without introduction of disturbance to the flow itself. 
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